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Abstract: At nonzero quark chemical potential dynamical lattice simulations of QCD 
are hindered by the sign problem caused by the complex fermion determinant. The severity 
of the sign problem can be assessed by the average phase of the fermion determinant. In 
an earlier paper we derived a formula for the microscopic limit of the average phase for 
general topology using chiral random matrix theory. In the current paper we present an 
alternative derivation of the same quantity, leading to a simpler expression which is also 
calculable for finite-sized matrices, away from the microscopic limit. We explicitly prove 
the equivalence of the old and new results in the microscopic limit. The results for finite- 
sized matrices illustrate the convergence towards the microscopic limit. We compare the 
analytical results with dynamical random matrix simulations, where various reweighting 
methods are used to circumvent the sign problem. We discuss the pros and cons of these 
reweighting methods. 
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1 Introduction 

In dynamical lattice simulations of QCD at nonzero quark chemical potential /i the gen- 
eration of a Markov chain through importance sampling is hindered by the sign problem 
caused by the complex fermion determinant, see [1] for a review. The severity of the sign 
problem grows as // increases and the determinant fluctuates more strongly. In the e-regime 
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of QCD, i.e., to leading order in the e-expansion of chiral perturbation theory [2], the spec- 
tral properties of the Dirac operator are universal and can be computed in the microscopic 
limit of chiral random matrix theory (chRMT) [3-5] . This equivalence also holds at 7^ 
so that chRMT can be used as a tool to investigate the sign problem. 

The fluctuating behavior of the fermion determinant can be characterized by is its 
average phase. Using chRMT, Splittorff and Verbaarschot have computed the average 
phase at / for trivial topology in the quenched and unquenched case [6]. Their 
results were later extended to nonzero temperature [7] and to general topology [8]. The 
complex analysis employed in ref. [8] is quite involved, and in the present work we give 
an alternative derivation of the formula for general topology, based on ideas presented in 
ref. [9] . Although the final integral expressions for the microscopic limit of the average phase 
look quite different in both cases, we show that they are indeed equivalent. En-passant 
we also derive some interesting new integral identities. In addition, the new derivation 
also provides an analytical expression for the average phase of finite-sized matrices, away 
from the microscopic limit. This allows us to verify the analytical formulas numerically 
using dynamical chRMT simulations. Such simulations are very costly and can only be 
performed with high statistical accuracy for small-sized matrices. In our dynamical chRMT 
simulations the complex weights are implemented using various reweighting methods. We 
provide a discussion of the pros and cons of these methods. 

The structure of this paper is as follows. In section 2 we introduce the chiral random 
matrix model with a chemical potential. In section 3 we show how the average phase of the 
fermion determinant can be computed in this model using complex Cauchy transforms. In 
section 4 the complex Cauchy transform is solved for finite-sized matrices, and in section 5 
the microscopic limit is taken. In section 6 we prove the equivalence of the integral repre- 
sentations of the microscopic limits derived here and in ref. [8]. In section 7 we verify the 
analytical predictions for the unquenched case by random matrix simulations away from 
the microscopic limit, using different reweighting methods. Finally we draw conclusions in 
section 8. Intermediate steps of the calculations are worked out in several appendices. 

2 Random matrix model 

Throughout this paper we use the same conventions as in ref. [8]. To make the presentation 
self-contained, we reproduce some of the equations derived in that paper. Details omitted 
here can be found in [8]. 

We work with the non-Hermitian chiral random matrix model for the Dirac operator 
D in the presence of a quark chemical potential introduced by Osborn [10], 



where the matrices ^pi and ip2 are complex random matrices of dimension {N + v) x N . 
They are distributed according to a Gaussian weight function given by 
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For a detailed analysis of this model, see also ref. [11]. Since the matrix in eq. (2.1) has \u\ 
exact zero modes we can identify u with the topological charge. From now on we assume 
u > 0; the results for v < follow by the replacement — t- in the final results. The 
nonzero eigenvalues of -D(^) come in N pairs {zj.,—Zk), and for /i = the Zk are purely 
imaginary. Note that /i in (2.1) is a dimensionless random matrix quantity and should not 
be confused with the physical chemical potential, see the beginning of section 5. 

For given i/, the partition function of the random matrix model with Nf dynamical 
quarks with masses is 

N f 

Zu^ {l-i;{mf]) = j dLpidip2w[vi)w{ip2)W_dei{D{ii) + nif) (2.3) 
with integration measure 

N+u N 

dX=lllldRe Xm d Im Xu • (2.4) 

k=i e=i 

The quenched case corresponds to Nf = 0. 

In ref. [10] it was shown that the partition function can be rewritten, up to a normal- 
ization constant, as an integral over the eigenvalues z^. of D, 



z!^^ {a; {nif}) 



N Nf 

Y[d^Zk w''{zk,zl;a) JJ("ij - zl) 
k=i f=i 



A^({.^})|^ (2.5) 



where we introduced q. = , the integrals over the z^ are over the entire complex plane, 

^n{{z^}) = W{zI-zI) (2.6) 

is a Vandermonde determinant, the weight function is given by 

..^(.,z*;a) = |zP-^exp ("^^^(-^ + ^*^)) (^^^1^1^) ' (2-7) 

and Ky is a modified Bessel function. The orthogonal polynomials corresponding to the 
weight function (2.7) are [10] 

Vl{z-a)={^-^^ k\Ll[-^^ , (2.8) 

where -^>^(^;) is the generalized Laguerre polynomial of order v and degree k. The corre- 
sponding orthogonality relation is 



Sz w'^{z, z*; a)p^{z; a)p£{z; a)* = r'^{a)5ke (2-9) 

c 



with norm 



Trail + a)^^+-kl{k + u)l 
= N^k+u+2 • (2-10) 
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The recurrence relation for the generahzed Laguerre polynomials, 

{k + l)L^+i(x) = (A; + 1 + u)LUx) - xL^^^x) , 
translates into a recurrence relation for the orthogonal polynomials p^, 

We will also use the Cauchy transform of the orthogonal polynomials defined by 



hUm;a) 



-,2 _ ^2 



w''{z,z*;a)p''^{z;ay 



The ensemble average of an observable O is given by 

- AT Nf 



1 



Y[ <fzk w''{zk,zl;a) Y{{mj - zl) 
k=i f=i 



(2.11) 



(2.12) 



(2.13) 



Ajv({^'})pO(zi,...,z^). 



(2.14) 



We will frequently omit one or both of the subscripts on (O). 



3 Average phase of the fermion determinant 

Adding a quark mass m to the Dirac operator we define D(m; fi) = D^jj) + ml, where m 
is assumed to be real. Writing detZ)(m;/u) = re*^, the phase of the determinant follows 
from [6] 



^2ie ^ det(j:>(/i) + m) 
° det(L>t(^) +m,) 



N 

n 



m 



(3.1) 



Here, m is viewed as a valence quark mass. We are interested in the ensemble average of 
e^*^ with two light dynamical quarks that have the same mass as the valence quark. This 
quantity is a measure of the fluctuations of the two-flavor determinant in the QCD weight 
function. For brevity we call e^*^ the phase of the determinant, although it really is the 
phase of the two- flavor determinant. Note that the average phase is real since each matrix 
appears in the ensemble with the same probability (2.2) as its Hermitian conjugate. 

In the presence of Nf dynamical quarks the average phase for a valence quark of mass 
m is given by 



^2ie\ _ I det{D{fi) + m) 
\det{D^{fi) + m) 



'Nf 



Nf 



Zf^(a;{m;}) ^c^ 

Zu' [a,m;{mf}) 
Z^f{a;{mf}) 



N 



Nf 



Jl (fzkw''{zk,zl; 



m 



a 



k=l 



2 



n(' 



\^n{{z'})\' 
(3.2) 
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where Z^^ is given by eq. (2.5) and Z"''^^^^ is the partition function of a random matrix 
model with iVj + 1 fermionic quarks and one conjugate bosonic quark, see ref. [6] for a 
detailed discussion. Both partition functions can be interpreted, up to a common additional 
normalization factor Z^, as averages of ratios of characteristic polynomials in the quenched 
ensemble. Such averages can be computed in terms of the orthogonal polynomials (2.8) 
and their Cauchy transforms (2.13) using the formalism developed in refs. [12, 13]. The 
details of its application to eq. (3.2) can be found in ref. [8, section 3.1], and one obtains 



PN-i{'m;oi) pj^_-^^{m;a) ■■■ p^_\ {m;a) 



PAr_l("lAf/;a) PAr_l("iAf/;a) ••• Vf^-\ \jnNj\OL) 



Nf 



U'^liimj - m2)J det [p"^^ \mf; a)] j g^^^ ^^ 
where we introduced the notation 



Nf 



(3.3) 



p/ {z;a) = z Pf^ {z;a) 



(3.4) 



and defined the complex integral 



w {z,z ■,a)p^_-^{z ) 



^•^_i(a) 7c z"^ - "T-^ 
over the orthogonal polynomials. In the quenched case (3.3) simplifies to 



•liB 



)nj=0 



p^_^{m;a) Pj:^_^(m;a, 



(3.5) 



(3.6) 



We now consider eq. (3.3) for the special case in which all dynamical fermions have the 
same mass m as the valence quark. We perform a Taylor expansion of the entries p'^'^{mf) 
of the determinant around m, 



p''i;^{mf;a) =p);^{m;a) + ^ — d^^p"/'' {m; a) , f = l,...,Nf. (3.7) 



A determinant remains unaltered when linear combinations of its rows are added to any of 
the rows. Therefore for each additional fermion it is sufficient to keep the next higher-order 
term in the expansion (3.7). The lower-order terms do not contribute to the determinant 
since they are identical to the contribution from one of the previous fermions, while the 
higher-order terms can be neglected since their contribution vanishes for my — t- m. After 
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taking each fermion mass in turn to m, we obtain 



dmP^N-iim; a) dmPN^im; a) 



■■ 'Hu,Nf+i{a,m) 
•• OmPN-i {m;a) 



Om PM_i[m;a) dm p j^j _i(jn; a) ••• dm p^-i ym;a) 



(2m)^/iV;! det [d!r^' p^^^-\m- «)] ^,,=1,...,^^, 
which can also be written in the form 

1 W^;\a,m) 



(3.8) 



{2m)^fNf \ Wl^^iO,l,...,Nf-l) 



(3.9) 



where 



Nf+l 



W^~\a,m) = {-)^n,,k{a,m)W^-_^\{0,...,k-l,k + l,...,Nf + l) 



(3.10) 



fc=0 



is a sum of Wronskians of order Nf + 1 with indices ranging from to Nf + 1, where in 
each term a different index k is absent. The Wronskian 



W^ip';'''{m- a), . . .yi'-{m; a)) = det [d'^^pf^im- a)] 
in eqs. (3.9) and (3.10) is abbreviated by W^{ki, . . . , kn)- 



(3.11) 



4 Solving the complex Cauchy transform 

The complex integral (3.5) needed in the computation of the phase factor is strongly os- 
cillating and cannot easily be evaluated numerically to high accuracy. In ref. [8] we solved 
this integral in the microscopic limit using quite involved complex analysis. We intro- 
duced integration contours which were then deformed such that the result was composed 
of contributions from the branch cut discontinuity and from the singularity of the modified 
X-Bessel function. The final result was given in terms of well-behaved one-dimensional 
integrals plus a short double sum. In this section we present a different derivation, which 
has the additional advantage of giving a calculable result for finite A^, away from the 
microscopic limit. 

We first explicitly substitute the weight factor (2.7) in eq. (3.5), 



'Hu,k{oi,m) 



1 



|^|2(.+l)^-a(.^+.-)^^(^|^|2)^^fc_^(^*.^)^ (4.1) 



where we defined a = N{1 — a)/4a and b = A^(l + a)/2a. In ref. [9], Osborn, Splittorff and 
Verbaarschot solved the Cauchy transform (2.13), which corresponds to the special case of 
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A; = in eq. (4.1).^ In the following we extend their result to arbitrary positive integer k, 
as required in eq. (3.3). Following the derivation in section V of ref. [9] we write 



+ (4-2) 



^2 _ ^2 -,2 _ j^2 ^2 _ j^2 

and decompose eq. (4.1) accordingly as Hiy^k = 'H^i^l + T~i^y \ with 

= / -^r^M'^''^'^^-^'''K^m^)p'^^_,{z*,a), (4.3) 



n^;fl{a,m) = -^^\- , % e— *V.(6|zp)p^^_,(z*;a). (4.4) 

For we consider the cases |z| < m and |z| > m separately and expand the deno- 
minator in a geometric series, 

' Af-lV"/ . n '''' J U <rri 



i=o 



oo „ 



j=0 



1 . 1 

(4.5) 



The polynomials p^''^(z;a) defined in (3.4) are even in their argument z so that we can 
implicitly define the expansion 

oo 

/(^2) ^ e— «) = ^ a„z2- . (4.6) 

n=0 

In polar coordinates the angular part of the integrands in eq. (4.5) can therefore be written 
as a product of power series in z and With z = re , the angular integration of such 
a product can be computed analytically using 

/ de z^^z*^^ = 2TTr*''6ki . (4.7) 
Jo 

After angular integration the second integral in eq. (4.5) vanishes, while the first one gets 
contributions 

/ de z'^^fiz*^) = V a„ / de z^^z*^"" = . (4.8) 

Jo Jo 

Resumming the aj using eq. (4.6), the integral (4.5) can then be written as 

duu''+^K^{bmu)e-'''''p''/_^{u; a) , (4.9) 



27re- 


am? 


/ dr 
Jo 




-i(a) 




-am? 


/ dw 
Jo 


' N-l 


(a) 



^We thank Jac Verbaarschot for drawing our attention to the solution of the Cauchy transform in ref. [9]. 
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where in the last step we have introduced the variable transformation u = r^/m and 
replaced /(u^) by its explicit expression (4.6). 

We now turn to n^ l. given by eq. (4.4). The second term on the RHS of eq. (4.2) is 
analytic and can be expanded as 



2 _„rr,2 OO 



-az „—am 



C 6 /I ^2 

n=0 



2-2 _ 



e— ^«(^;«), (4.10) 



where the coefficients dn were computed in ref. [9] and are given by 

dn = - „ , / dte-'"'^*/^!-"), • (4-11) 

Note that the explicit factor e~"^ in eq. (4.10) is such that the weight of the orthogonal 
polynomials are retrieved after substituting eq. (4.10) in eq. (4.4). This yields 

'H^^l{a,m) = — i—,^dn I cfzw''{z,z*;a)p';^{z;a)p''^ti{z*;a). (4.12) 



In ref. [9] the integral (4.12) was solved for k = 0. In that case the solution immediately 
follows from the orthogonality relation (2.9), resulting in 

, , ]\rN /.(l-a)2/4a 

nl^Ua,m) = -d^., = / rf,,->n-A^V(i-^) , (4.13) 



(A^-l)!(l-a)^ Jo (t + l)^+'^ 

where we substituted the explicit expression for d]^^i given in eq. (4.11). To compute 
the phase of the fermion determinant for arbitrary Nf and z/ we also need the solution of 
eq. (4.12) for positive integer k. In this case the second polynomial in the weighted integral 
is p^j^_^{z*; a), which involves an orthogonal polynomial of order v + k instead of v and an 
additional power of z*, such that the orthogonality relation (2.9) can no longer be applied 
directly. To solve the integral (4.12) for arbitrary /c G N we use the relation 

pfiz-^a) = z^-f-\z,a) = (4-14) 

for any v,^,k G N, which expresses the LHS as a sum over orthogonal polynomials of 
order v and degree £,...,£ + A;. The proof of this relation is given in appendix A. We now 
substitute this expansion in eq. (4.12) to find 

k 



(2) ^ 1 (k\ {N-l + u + k)\ (a-lV 

"^'^^"'""^ r-^_M)jT;,\j)iN-l + v + k-j)\ \ N ) 

oo „ 

y^dn d^zw''{z,z*;a)p'',Xz]a)p''j^^i^,,_Az*]a) 

Jc 

h U {N-i + u + k-jy. [-it) ' 

(4.15) 



X 

n=0 
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where in the last step we have apphed the orthogonahty relation (2.9), after which only 
the n = N— 1 + k — j term survives. Using the norm (2.10) we compute the ratio 



rXr-i+fc-,-(«) ^ (1 + a)2(fe-J-)(jV - I + k - j)\{N - 1 + k - j + 1^)1 

r5^_i(Q) iV2(fc-i)(Ar_i)!(Ar_i + j,)! ' ^ ' ^ 

and after substituting this and the expression (4.11) for d in eq. (4.15) we find 



(1 - 0)^+^= {N - 1)\{N + Jo {t + l)N+k+u 

{N-l + v + k)\ 



3=0 



{N - iy.{N - 1 + i^y. 



1 . n^-ifu-l^'^ 



due-"" """j- ^ , . (4.17) 

[(1 - a)2n + 4a]^+'=+'^ ^ ' 

In the last step we have performed the sum over j using the binomial theorem to obtain 
[Aat — (1 — a)"^]^. Also, we have introduced the variable transformation u = 4at/(l — a)"^ 
to simplify the integration range. 

The sum of H^j^l and T-L^i^l of eqs. (4.9) and (4.17) gives the exact finite- result for 
T-Lu,k needed to compute the average phase of the fermion determinant with eqs. (3.3), (3.6) 
and (3.9). For < one just needs to replace by Note that the chiral limit (m — t- 0) 
of the phase factor (3.9) has to be taken carefully, as detailed in appendix B.l. 



5 Microscopic limit 

The results computed from chRMT are universal, i.e., identical to the corresponding quan- 
tities in QCD, in the so-called microscopic regime. This regime is obtained by defining 
the rescaled parameters rh = 2Nm, rhf = 2Nmf, a = 2Na and the rescaled eigenval- 
ues z = 2Nz and then taking N ^ oo while keeping the rescaled quantities fixed. The 
conversion of the rescaled random matrix parameters to the physical parameters is done 
using the relations z = ^phys 

VT., rh = niphysVT. and a = fi'^ = ^'^^^^F'^V, where V is 
the four-volume and S and F are low-energy constants of chiral perturbation theory. The 
Gell-Mann-Oakes-Renner relation m'^.F'^ = 2mS (for equal quark masses) can be used to 
introduce the physical pion mass 771,7^ through the combination fip^^^/rri^ = fj?/2rh. 

To take the microscopic limit of the average phase we introduce the corresponding 
limits of the various objects defined in section 2 and section 3. A detailed derivation of 
these limits can be found in ref. [8, appendix A]. The microscopic limit (denoted by a 
subscript s) of the orthogonal polynomials (2.8) is defined as 

p^(i; a) = lirn^ —f—^p-^_^(z/2N- a/2N) = ^e-^'H-^I, (z) . (5.1) 
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Accordingly, the microscopic limit of p^]^_i, defined in eq. (3.4), is 



p^s^{z\ a) = lim 



pf_,{z/2N;a/2N) = ^e-^'^m^^k {z) 



where we introduced the notation 

Iv,k{z) = z''I^+k{z) . 

The microscopic limit of the weight function (2.7) is 

<(z,r;a) = lim {2Nf''+^w^{z/2N, z* /2N;a/2N) = |i|2('^+i)e-%^K^ , 

7V-s>oo y 4q! 

and the microscopic limit of the normalization factor (2.10) is 



(5.2) 



(5.3) 



(5.4) 



(5.5) 



<(a) = lim (2A^)"e"^V^„i(d/2A^) = Air'ae" . 

We define the microscopic limit of the integral (3.5) as 

ntk(.a,m)= lim e-^{2Ny+''-'^/^V^e-^/^m-'"H^k{a/2N,rh/2N), (5.6) 

where the prefactors were chosen such that it coincides with the master integral of ref. [8, 
eq. (3.21)] and remains finite when N — t- oo, i.e., 



-2a 



9 o -e Kt, { ' — -]Ii,kiz*)- 

Anarfi'^ r ^ — z*^ V 4a ' 



(5.7) 



Here we also used the microscopic limits (5.2), (5.4) and (5.5). Using eqs. (5.2) and (5.6) 
we now construct the microscopic limit of the phase factor derived in section 3. For the 
unquenched case we take the microscopic limit of eq. (3.3) and find 

Iufi{m) Iu,i{m) ■■■ Ii,^Nf+i{rn) 
Iufi{rni) Iu,i{mi) ■■■ Iu,Nf+i{'mi) 



2ie 



lim {e''')^ 



XffUim) - m2)J det [/^,g_i(m/)] 



(5.8) 



As expected, the dependence on has dropped out, leaving a finite microscopic limit for 
the average phase factor. Similarly, the microscopic limit of the quenched average phase 
factor (3.6) is given by 



hm (e^^^)^ _„ 



Iufl{m) Iy,i{rh) 



(5.9) 
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For the special case in which all dynamical fermions have the same mass m as the 
valence quark, eq. (5.8) becomes 



/Nf 



Iu,o{in) Iu,i{rh) 



I, 



u,Nf+l 



lire-) 'ir 



m] 



(m) 



i2m)^fNf\ det [liiZ'H^)] f,g=i,...,N, 
in analogy to section 3. Again, an alternative way to write this result 

'^f {2m)^fNf\ 1, . . . , TV/ - 1) ' 



(5.10) 



IS 



(5.11) 



where we have defined 

Nf+l 



WNf{a,m)= ^ (-)'=?^^^fc(Q,m)VF^^+i(0,...,A;-l,fe + l,...,A^/ + l) (5.12) 

k=0 



in analogy to (3.10) and Wn{ki, . . . ,kn 



is a short-hand notation for the Wronskian 
WMuMi^), ■ ■ ■,Iu,kA^)) = det [4X'^(rn)] (5.13) 

appearing in eqs. (5.11) and (5.12). 

To compute the microscopic limit (5.8) of the average phase we need to take the 
microscopic limit (5.6) of the complex integral (3.5), which was solved for finite N by 
eqs. (4.9) and (4.17). Taking the microscopic limit (5.6) of eqs. (4.9) and (4.17) gives 



n - "1 



4c. 7. ^-'=-^^^-«"^^^(^)w.(n) 



+ \, : .. I due u 8& ^ ^ 



where in the first term we have substituted eqs. (5.2) and (5.5) and in the second 
have used the definition 



lim 1-^ = 



Af^oo \ N J 

of the exponential function and Stirling's formula 



= e 



-a 



(5.14) 
. term we 

(5.15) 



m 

lim , = 1 . 



For the special case of A; = 0, eq. (5.14) reproduces the result which was computed pre^ 
viously in ref. [9] in a study of the chiral condensate. ^ — — „ 



(5.16) 
pre- 

However, for the calculation of the 
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average phase using eq. (5.8) we need the more general result of eq. (5.14) for arbitrary 
k gN. The chiral limit of the phase factor (5.11) is computed in appendix B.2. 

In ref. [8] we computed the complex integral (5.7) in a different way, which led to an 
integral representation that looks quite different from eq. (5.14). In the next section we 
will prove the equivalence of the expressions derived in both formulations. 

Apart from providing an alternative formula to compute the microscopic limit of the 
average phase, the current derivation has the significant additional feature that it also 
gives a calculable and well-behaved expression for finite-sized matrices, away from the 
microscopic limit, given by the sum of eqs. (4.9) and (4.17). This expression will be useful 
when verifying the analytical RMT results by dynamical RMT simulations in section 7. 



6 Equivalence of integration results 

In ref. [8] the complex integral (5.7) was solved using a completely different formalism, 
based on the deformation of integration contours and some involved complex analysis, 
leading to 



— 2a— 

4a 



mu 
4a 



Iu+k{u) 



3.1) 



+ 



(iu(-)V+ie-8S/^ 



mu 
4a 



Ku+kiu) 



+ e 8a ^-L— 

rh 



— E 



{u-l-i)\{v-l + k-j)\ (nl_V J 



This solution must be equivalent to (5.14), as they are just two different representations of 
the same complex integral (5.7). Nevertheless, as both expressions look quite different it 
is useful to prove their equivalence, without resorting to the complicated complex analysis 
used in ref. [8]. This proof will also provide a check on both results. 

The first integral in eq. (5.14) is identical to the first integral in eq. (6.1). Therefore it 
remains to show that the second integral in eq. (5.14) equals the sum of the second integral 
and the additional polynomial in eq. (6.1). Setting a = 2a and b = rn^/Sa to simplify the 
notation and canceling some prefactors, we thus need to show 



-bu 



-a—b 



W 



v+k+l 



' duv!'''^e'^Iu{u^/h/a)K^+k{u) + Su,k 



2^a 



3.2) 



where we have defined 

Sv,k 



(i/-l-^)!(z.-l + fc-j)! ^,^,- 



i,i=0 



(v-l 



j)!i!j! 



(6.3) 



For the special case oi v = k = {) this was proven in ref. [9]. Below we will give a general 
proof for arbitrary v and k. 
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In appendix C we prove the integral relation 



OO ,,2 

k+1 



duu ~^ e ^aIy[u\Jh/a)Ky+k{'u) 



u/2 



ds e 



(6.4) 



Renaming u to s on the LHS of eq. (6.2) and substituting (6.4) in the first term on the 
RHS, proving (6.2) corresponds to proving the identity 



ds e 



gU+k+l 



6 Sy k ■ 



Using the Rodrigues formula for the generalized Laguerre polynomials, 

e^x'" d^ 



Ll{x) 

the LHS of eq. (6.5) can be written as 



X,fe = (-)V / dse 



k\ dx^ 



(6.5) 



(6.6) 



(1 - s)" ds'' 



[e-'^l-s) 



ds 



ds'^ 



e s 



[i-sY 



e-'^l-s) 



(6.7) 



where in the last step we have performed k successive integrations by parts. In appendix D 
we prove that 



d^ 
ds'' 



e s 



{i-sY 
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e~ s 



j=0 



k\ {y)ja ^s 



k-j u-k+j-1 



jj (1 - s)^+i 



with the Pochhammer symbol defined in appendix C, and substituting this in the 
integral (6.7) we obtain 



(6.9) 



With this identity we show in appendix E that Ju,k satisfies the recurrence relation 



X,fc = -b''-^e-''-\iy)k + ^{j + l)k^j [a6X-2,i + (i^ - l)X-ij] • (6.10) 

j=0 

We now prove the equivalence (6.5) by induction in using the recurrence relation 
(6.10) and the expression (6.9) for Ju,k- For 1^ = only the j = term contributes in 
eq. (6.9), and we find 



Jo,k = a'' I dse~s~''s~'"\l-s)\ 



(6.11) 
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which corresponds to the RHS of eq. (6.5) since S^fc = 0. For v = 1 eq. (6.9) gives 
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k-j 



a k-j 



s(l-s) 



kl 



.k-j+ 



,+1 (1 - sf-^ 



gk-j+2 



t e s 



„fc+2 



(6.12) 



where in the second step we have integrated by parts and in the third step we have observed 
that the j = k contribution in the last term vanishes and shifted j ^ j — I for this term 
so that in the last step only the j = term of the first sum remains. For u = 1 eq. (6.3) 
gives Si^k = k\ so that eq. (6.5) is satisfied for ly = 1. We now prove that the equivalence 
holds for J^^k if it is satisfied for J'u~i,k and J'u~2,k- We start from the recurrence relation 
(6.10) and substitute the RHS of eq. (6.5) for — 1 and 1^ — 2, 



k 

X,fc = Va-~i+^(i + l)fc_, / 
,=0 -^0 



1 a , (I — o)j 



(6.13) 



Performing an integration by parts on the first term in the integral yields 
-1 







ds [ 



" " — ' dse s 







1 — s s 



„ — a—bj; 



(6.14) 



where we observed that the surface term is only nonzero for j = 0. Let us denote the first 
term on the RHS of eq. (6.13) by JT^^''. Substituting (6.14) into this term we obtain 



j=0 

k k-1' 

E-E 

b=o j=o. 



\ -^-bs {i-sy 



s"^ s(l — s) 



e-^-^a^-^kl 



f 

JO 



7 — — —OS 

dse s 



^v+k 



(6.15) 



In the first line, we observed that the j = contribution to the second term in square 
brackets vanishes, which allowed us to shift j — >• j ' + 1 in the corresponding sum in the 
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second line. In the third line only the j = k term of the first sum survives. Putting 
everything together, we have 



Jo 



(l-g)fe 

gu+k+l 



- + [abS,.2,j + l)Su-i,j] + b'-Hu);, + a'^-ifclj . (6.16) 

To reproduce the RHS of eq. (6.5) it remains to show that 
k 

J2U + ^)k-j [abS,-2,j + 1)5.-1,,] + b''-H^)k + a'-^kl = S,,k , (6.17) 
i=o 

which we relegate to appendix F. This completes the proof of the equivalence (6.5) for 
arbitrary i' and k. 



7 Dynamical random matrix simulations with complex weights 
7.1 Computing averages using reweighting 

In ref. [8] the analytical chRMT results for the microscopic limit of the quenched phase 
factor were thoroughly verified using quenched random matrix simulations. In this section 
we verify the unquenched analytical results using dynamical random matrix simulations. 
In the dynamical case the finite- corrections to the microscopic results are quite signif- 
icant. Increasing the matrix size in the simulations to be close to the microscopic limit 
would require too much computational power, especially since dynamical simulations are 
already intrinsically expensive. Instead, we chose to use small-sized matrices and compare 
the results for the average phase factor with the finite- predictions of eq. (3.9), where the 
master integral is given by the sum of eqs. (4.9) and (4.17). Even then, the random ma- 
trix simulations at nonzero chemical potential are problematic since the partition function 
contains a complex weight. Both the real and imaginary parts of the weight function can 
become negative, and as a consequence the fermion determinant cannot be included in the 
probability distribution of a Markov Chain Monte Carlo (MCMC) simulation, i.e., we are 
confronted with the sign problem. 

To circumvent this problem we will perform the importance sampling using an auxil- 
iary, nonnegative, weight function, and reweight the result appropriately so that we average 
over the correct target ensemble. This reweighting leads to an overlap problem, when the 
configurations contributing most to the partition functions in both ensembles do not co- 
incide. Note that this overlap problem always occurs in reweighting, even if the weight in 
the target ensemble is nonnegative. This is, for example, also an issue at zero chemical 
potential when using dynamical fermions with a heavier mass in the auxiliary ensemble 
than in the target ensemble. In the presence of a chemical potential, the overlap problem 
is amplified by the sign problem. Even if the auxiliary distribution has a good overlap with 
the target ensemble, i.e., the relevant configurations are sampled appropriately, the sign 
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problem can ruin the reweighting procedure. This wih happen due to large cancellations 
of contributions with opposite sign in the partition function of the target ensemble, which 
will generate large statistical errors. 

Below we briefly describe the principle of reweighting.^ Using the partition function 
(2.3) we want to compute the unquenched expectation value 



{0)Nf = ^;[j j d<^idLp2w{Lpi)w{Lp2)'D[ipi,Lp2]li\{mf]]0[Lpi,^p2] (7.1) 
with dynamical determinant 

V[^i, ip2; fi; {rrif}] = J] det(D(^) + mf) = Re'® . (7.2) 
/=i 

Here, R is the (nonnegative) magnitude and e*® is the phase of V. In the case of interest 
for this study, the observable will be O = e^*^. To set up the reweighting formalism we 
introduce the weighted (or ensemble) average 

_ Jdx w{x)f{x) 

\J/W — r , / N I'-Jj 

J ax w[x) 

of / with respect to w, with normalization (1)^ = 1. From this definition we see that 
an ensemble average can be computed from an auxiliary ensemble using the reweighting 
relation 

{gUf = ^ . (7.4) 

This feature is useful to study ensembles with weight functions that cannot be sampled 
efficiently, or where the weights are not positive definite such that they cannot be used 
as probability distributions in importance sampling. The actual simulation constructs a 
Markov chain for an auxiliary ensemble, after which the expectation value of the observable 
in the target ensemble is computed by reweighting the observable and the partition function 
as given in eq. (7.4). To keep the statistical error of the reweighted observable within 
reasonable limits, the overlap between both ensembles should be large, i.e., the bulk of 
relevant configurations in both ensembles should coincide. 

In this study we compare the results obtained with three different reweighting schemes 
[19]. The ensembles and the corresponding reweightings are as follows. 

Rl Quenched simulation with full reweighting: We perform a standard quenched simula- 
tion of random matrices through direct sampling of the Gaussian weights (2.2) for the 
real and imaginary parts of the elements of ipi and cp2, as described in appendix E of 
ref. [8], and reweight with the dynamical fermion determinant (7.2). Using eq. (7.4) 
the unquenched average is computed from 



^For related approaches see, e.g., [14-18]. 
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where the quenched averages in the numerator and denominator are evaluated as 
averages over the Monte Carlo sample. 

The omission of the fermion determinant in the importance sampling will cause an 
overlap problem, and the reweighting factors will strongly fluctuate between config- 
urations. These fluctuations will increase as the chemical potential grows and the 
sign problem becomes more tangible. In this reweighting scheme the generation of 
the matrices in the Markov chain is cheap, but we need a very long chain to get an 
acceptable accuracy. 

R2 Phase quenched simulation using Metropolis with partial reweighting: In this Metro- 
polis algorithm the matrix probability distribution consists of the product of the 
Gaussian weights (2.2) and the magnitude R of the dynamical determinant. The 
measurement is then reweighted by the phase factor of the dynamical determinant, 

Including information about the determinant in the sampling probability should im- 
prove the overlap between the generated configurations and the significant configu- 
rations in the unquenched ensemble. 

R3 Sign quenched simulation using Metropolis with minimal reweighting: As the un- 
quenched partition function (2.3) is real, the contributions from the determinant 
(7.2) to the partition function only come from ii cos 0. The imaginary contribu- 
tions cancel between Hermitian conjugate matrices as these have the same Gaussian 
probability (2.2). To minimize the variance of the reweighting factors and optimize 
the overlap it is therefore natural to create an auxiliary ensemble using the weight 
R\ cos Q\ and absorb the remaining sign of the determinant in the reweighting factor, 

un\ ((sgn cos 6) (1 + i tan 9) 0)^1^08 01 „s 
^ (sgncose)/j|cose| 

where the denominator is explicitly real because of the symmetry mentioned above. 
Using the absolute value of the weights as auxiliary distribution allows one to include 
as much information as possible about the determinant in the MCMC weights. Note 
that the average over the itanQ term in the numerator does not vanish in this 
case because the observable O = e^*^ also has an imaginary component. We expect 
this reweighting scheme to be somewhat more effective for real than for complex 
observables. 

After performing simulations using this reweighting, we realized^ that this idea of 
minimal reweighting had been discussed earlier [20] in a study of the reweighting 
factor in lattice QCD simulations at nonzero fi. In that study is was shown, based 
on the central limit theorem, that reweighting by i?|cos0| indeed minimizes the 
fluctuations in the reweighting factor. In ref. [20] this reweighting scheme was not 



^We thank Philippe de Forcrand for bringing this earher study to our attention. 
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implemented since it would have been too expensive in lattice QCD. Here, we put 
it to a practical test on the random matrix model. (Note that scheme R3 was also 
rediscovered in [21].) 

Note that for a complex observable one can develop more sophisticated reweighting 
algorithms using two Markov chains for the numerator, one for the real part of the 
weights and another for the imaginary part. The complication is that the chains have 
to be normalized properly with respect to each other, which introduces additional 
overhead. This will not be investigated further here. 

Another partial reweighting scheme uses the //-quenched ensemble, where the matrices 
are sampled using a dynamical MCMC algorithm at zero chemical potential, and the 
contributions to the averages are reweighted by the ratio of the determinants at chemical 
potential // and at /i = 0. For small /i, this scheme is close to scheme R2, and for larger fj, 
it is no better than scheme Rl but much more expensive. Therefore it will not be studied 
here. 

We briefly discuss the main reweighting features of the three schemes. Clearly, the 
overlap problem is best handled by the sign quenched scheme, which will sample the most 
significant configurations of the partition function. The biggest overlap problem will be 
encountered by the quenched scheme, as the fermion determinant is completely ignored in 
the auxiliary weight function. Nevertheless, the construction of the Markov chain is very 
cheap in the quenched scheme, so that the overlap problem can be partially alleviated by 
the generation of many more configurations, which are all uncorrelated by construction. 
Both the sign quenched and phase quenched schemes use a Metropolis algorithm. Here 
the autocorrelation time has to be taken into account to select the independent configura- 
tions, which very much shortens the effective size of the Markov chain. Fortunately, the 
much smaller fluctuations allow us to reach a high accuracy with much fewer uncorrelated 
configurations than in the quenched scheme. Although the sign quenched scheme seems 
superior to the phase quenched one in the case of random matrix simulations, the latter can 
be more easily implemented in realistic theories like QCD. A phase quenched determinant 
for Nf = 2 can be implemented by simulating a quark and a conjugate quark, whereas 
a sign quenched determinant does not seem to have a physical equivalent which could be 
implemented efficiently. 

Even though the severity of the overlap problem is different for the three schemes, 
there is no reason why the sign problem would be improved upon in any of the schemes. 
Indeed, even in the sign quenched case where we sample the most significant configurations, 
the positive and negative contributions will balance each other more and more when the 
chemical potential becomes large, and the sign problem will remain. This will be confirmed 
in the numerical experiments discussed in the next section. 

In the numerical implementation, for each matrix generated in the importance sampling 
we also consider its Hermitian conjugate matrix, in accordance with the symmetries of the 
partition function. This ensures that the sample average of the phase factor is explicitly 
real, which somewhat simplifies the implementation of the algorithm and the computation 
of the statistical error on the final result. 
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Figure 1. Average phase factor of the fermion determinant for Nf = 2 as a function of the 
fermion mass rh for a — 1.0 and v = 1 for matrix sizes A'^ = 2, 4, and 8. The simulation results, 
computed with sign quenched reweighting, are given by the data points, while the solid lines show 
the corresponding analytical results of eqs. (4.9) and (4.17). We also show the microscopic limit 
(N ^oo) of eq. (5.14). 

7.2 Numerical results 

We performed dynamical random matrix simulations at nonzero chemical potential using 
the three reweighting schemes described in the previous section. For each measurement 
we generated 1,000,000 matrices with sizes ranging from = 2 to 16. In the quenched 
ensemble the generated matrices are uncorrelated, but for the sign and phased quenched 
ensembles, which are sampled using a Metropolis algorithm, successive matrices in the 
Markov chain are correlated. This autocorrelation effectively reduces the number of in- 
dependent configurations and is taken into account appropriately when computing the 
statistical errors on the measurements. In addition, the statistical errors take into account 
the correlations between numerator and denominator in eq. (7.4). Details of the calculation 
of the statistical errors are given in appendix G. 

In figure 1 we verify the mass dependence of the phase factor for N = 2, 4, and 8 for 
Nf = 2 and v = 1. The simulation results agree very well with the analytical predictions 
of eq. (3.9). Both the m-dependence and the dependence on the matrix size A^ of the 
matrices is reproduced. The different curves show how the large-A'" limit is approached. 
For larger masses the microscopic limit is not yet reached for matrix sizes up to A^ = 8, 
and the analytical results for finite A^, away from the microscopic limit, are essential to 
explain the numerical results. 

Figure 2 illustrates how the phase factor changes as a function of the chemical potential 
for matrix sizes A^ = 2, 4, 8, and 16. The results again agree very well with the analytical 
predictions. To probe the overlap and sign problems, which manifest themselves in the 
accuracy of the measurements, the simulations were performed using the three reweighting 
algorithms discussed above. For small chemical potential the accuracy of all three algo- 
rithms is satisfying, reflecting the fact that neither the overlap nor the sign problem is 
significant. However, the phase and sign quenched algorithms achieve this accuracy with 
much fewer independent configurations than the quenched algorithm, as a consequence of 
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Figure 2. Average phase factor of the fermion determinant for Nf = 2 as a function of the 
chemical potential a for rh — 1.0 and ^ = 1 for matrix sizes = 2, 4, 8, and 16 using three 
different reweighting schemes. The blue squares were computed from a quenched ensemble with 
full reweighting, the red circles from a phase quenched simulation with partial reweighting, and the 
purple triangles from a sign quenched simulation with minimal reweighting. 

the better overlap between target and auxiliary ensemble. The shorter trajectories would 
be an advantage if the measurement itself were expensive. With increasing chemical poten- 
tial the error on the measurements increases as the sign problem sets in (for <^e^*^) < 0.2). 
The sign problem equally affects the different reweighting methods. As expected, the sign 
problem gets worse for larger N . In the Metropolis algorithm this deterioration for in- 
creasing N is related to an increasing autocorrelation time. Even for the sign quenched 
algorithm, where the overlap problem is minimized, the sign problem remains as the matri- 
ces are reweighted by +1 or —1, which leads to important cancellations and large statistical 
errors for large a. 

8 Conclusions 

Dynamical lattice simulations of QCD at nonzero baryon density are hindered by the sign 
problem caused by the complex fermion determinant. To investigate this problem it is 
helpful to employ the equivalence between the spectral properties of the Dirac operator 
in the e-regime of QCD and chiral random matrix theory, which also holds at nonzero 
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chemical potential. As the average phase factor of the fermion determinant is an important 
observable in the study of the sign problem, we have computed it in the framework of chiral 
random matrix theory. 

In ref. [8] we derived an analytical formula for the average phase factor in the micro- 
scopic limit of quenched and unquenched chiral random matrix theory for general topology. 
In the current paper we presented an alternative derivation, leading to a different integral 
representation of this microscopic limit, and showed that both formulations are equivalent. 
In contrast to our previous work, the new formulation also gives calculable expressions for 
finite-sized matrices, away from the microscopic limit. 

The analytical predictions of the finite-A^ formula were verified in dynamical RMT 
simulations, where reweighting techniques were used to compute averages in an ensemble 
with a complex weight. Very good agreement was found. At large chemical potential the 
statistical errors grow, signifying the emergence of the sign problem. 
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A Relation between orthogonal polynomials 

In this appendix we prove the identity (4.14) by induction. For A; = (4.14) trivially holds 
for any i^, £ S N. For integer A; > the recurrence relation (2.12) yields 



where we have first shifted the index j by one in the second sum and then gathered the 
overlapping terms in both sums. As the binomial coefficients satisfy 





Assuming that (4.14) holds for p'^'^ ^ and any iy,£ £ N, we can substitute it twice in the 
previous equation, for p^^^ ^ and p'^'^ ^ , to find 
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(A.3) 
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we find 

k 



j=0 

which shows that the identity (4.14) holds for arbitrary integer A; > 0. 



B Chiral limit 



B.l Finite N 



The chiral limit of the average phase (3.9) has to be taken carefully such that the mass 
factors in numerator and denominator are canceled properly. For this we need to compute 
the limit m — t- of the Wronskian (3.11), which contains derivatives of the function p'^'^ 
defined in eq. (3.4). For small argument m the leading-order term of this polynomial is 

^fc fl-aY {e + u + k)\ 2k 

f ™ 

with p-th derivative 



. 71- a^ (''+'' + *)' (2*)! 2k-p 



Substituting these expressions in the Wronskian (3.11) and using properties of the deter- 
minant gives the leading-order result 

(B.3) 

where A„(/ci, . . . ,kn) is a Vandermonde determinant. From this we compute the chiral 
limit of the denominator in eq. (3.9), 

<^(0,l,...,iV,-l)~ (^) (2m)^/(^/-i)/2 n ^ . +y n 
where we used the identities 

Nf~l Nf-1 

AM,.{0,l,...,Nf-l)= II £\ and i = -Nf{Nf - 1) . (B.5) 

1=1 i=0 

From eq. (B.3) it is easy to see that in the limit m — )• only the Wronskian corresponding 
to the k = Nf + 1 term in eq. (3.10) contributes to leading order, while all other terms are 
of higher order. The average phase (3.9) can therefore be written as 



\ /m=o m^o(2m)^fNf\Wj^^{0,l,...,Nf-l) ^ 
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where in the last step we have substituted (B.4). It now remains to take m = in ?^ 
computed using the integrals (4.9) and (4.17). The first integral is zero as the integration 
range vanishes, and hence the average phase simplifies to 



X'^ /m=0 
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{Nf + 1)\{N + Nf + uy. 
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(B.7) 



where in the second step we have used the integral representation of the hypergeomet- 
ric function [22, eq. (15.3.1)] and in the last step we have applied the transformation 
2Fi(a, b; c; z) = (1 - z)^-'^-*2i^i(c -a,c-b;c; z) [22, eq. (15.3.3)]. For > 1 the hyperge- 
ometric function is a polynomial of degree v — 1, see eq. (B.13) below, and for v = \ the 
result further simplifies to the A^j-independent expression 

2N 



m=0 



1 



a 



1 + a 



B.2 Microscopic limit 

We now take the N ^ 00 limit of the results in appendix B.l with d = 2A^a fixed. In this 
limit we can replace (1 + a)^ in the second factor and (1 — a)^ in the fourth factor of (B.7) 
by 1. We also have 

2N 



1 



a 



lim 

V 1 + a 



-2a 



(B.9) 



For = we use the relation 



N 



lim 2F1 {l,n:N- = lim (N - 1) 

7V-5>oo V X J N^oo 



dt- 



[1-t) 



N-2 



1 + Nt/xY 



dz- 



= 3;"e^r(-n + l,2;), (B.IO) 

where in the first step we have used the integral representation of 2-^1 [22, eq. (15.3.1)], in 
the second step we have substituted z = x + Nt and taken the N ^ 00 limit, and in the 
third step we have used the integral representation of the incomplete gamma function [22, 
eq. (6.5.3)]. With (B.IO) we obtain from (B.7) 

- [ ^ + ' ' 

_N^oo 2a N + Nf + 1_ 



X'^s /u=0,m=0 



{2a)^f+^e^^T{-Nf - l,2a) 



{Nf + l)(2d)^^+^r(-A^/ - 1, 2d) 



(B.ll) 
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For 1/ = 1 we obtain immediately from (B.8) 

(ef ) ^ , , = e-2". (B.12) 

For 1/ > 2 the first argument of the hypergeometric function in (B.7) is a negative integer, 
in which case we have [22, eq. (15.4.1)] 

2Fi{-m,b;c;z) = — (B.13) 

n=o 

and obtain 



/^2^e^ ^ -2a ( iV/ + 1) ! ( iV + iVj + z.) ! ^ (-r(l-i^)„(iV; + 2)„. f2a\ 
/'^=o N-^^{Nf + i^)l{N + Nf + l)\^^ nl {N + Nf + 2)n \nJ 



^ ' (A^/ + i/)!n!(z^-l-n)!^'^"^ nToo {N + Nf + n + ly.N"-'^- 



n=0 
v-l 



where in the last line we have substituted j = v — 1 — n. 

Eqs. (B.ll), (B.12) and (B.14) agree with the results previously computed in ref. [8]. 

C Proof of integral relation (6.4) 

Below we prove the integral relation (6.4). In the following, the LHS of eq. (6.4) will be 
denoted by Ty^^, i.e.. 



Iy,k= / duu''+^e-^Iy{u./b/a)Ky+k{u) . (C.l) 
Jo 



Using the integral representation [23, eq. (8.432.6)] 



^^(") = ^(i) / ^^^^ ifRe(n2)>0, (C.2) 

-' 



which is vahd for arbitrary z^, we obtain 



-| roo f'OO t-\-a 

The integral over u is given in [23, eq. (6.631.1)] in terms of a confluent hypergeometric 
function, 

T(u + k + l) /bV^^ / 4at Y+''+^ ^( , bt 



2-+ir(.+i) y-aj yiT-aj ^^-^^ 
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Using the Kummer transformation [22, eq. (13.1.27)] 

iFi(a; b; z) = e\Fi{b - o; 6; -z) 

we rewrite (C.3) as 



(C.5) 
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r(^^ + i) KaJ Jo 



(C.6) 



where in the last step we have introduced the variable transformation s = a/{t + a). 
Note that iFi {—k; + 1; z) = k\L^[z)/{v + 1)^ with the Pochhammer symbol {a)n = 
a{a + 1) . . . (a + n — 1) and (a)o = 1 so that the integral finally becomes 



I^^k = 2''k\a[ y e"+^y dse 



,-+'=-iL-(-6(l-.)). 



This proves the integral relation (6.4). 
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D Proof of identity (6.8) 

In this appendix we prove the identity 
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by induction. For A; = the identity is trivially satisfied for any i/, with 



For A; > 1 we perform one derivative in D^, ^ explicitly and obtain 



D 



ds^-^ 

d^ 

ds^-^ 



ae s ^ — \- (h' + k — l)e s ^ — \- i^e s ^ — — 



(D.l) 



(D.2) 



ae 



{i-sY 



^ + {a + v + k-l)e 



{i-sY 



+ ve s 



{l-s) 



v+l 



= aDi,^i^k-i + {a + ly + k- l)Di,^k-i + I'Du+i^k-i 
Assuming that eq. (D.l) holds for fc — 1 we thus have 

(i/ - l)ja^-^s''-'=+^-i 



(D.3) 



k-i ^ 

a /k — 1 



u,k 6 



E 

j=0 



J 



+ {a + u + k-l) 



k-j-l oi^-k+j 



(1 - s)^+J 



+ 



(z^)j+ia^-J-is^-'=+^'+^ 
(T^^p+j+i 
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e s 



■k-i 

E 

fc-i 

+ E 

i=o 



A; - 1 



(1 - sY+i-^ ^ (1 - sY+i 



k-1 



-,k-j-l 



+ 



(1 - sY+^ (1 - sY^^+^ 



a I ^ / A; - 1 



+E 



(Z^ - l)jS 



u—k+j 



k- 1 



(1 - sY+^-^ (1 - s)''+J' 



^ (1 - sY+^~^ (1 - s)^+J 



(1-.) 



/c — 1\ ( {v — l)jS 



J 



+ 



u-k+j 



(1 - s)'^+J-i (1 - 



+ 



(1-.) 



v+k 



(D.4) 



where in the second step we have gathered the terms with equal powers of a, in the third 
step we have shifted the index of the second sum, and in the last step we have extracted 
the j = and j = k terms. The binomial coefficients satisfy 

k /k-V 



k-l\ fk-l 



) and ( 

jj \jj j \j - 1 



and the Pochhammer symbols obey the identity 

u-l+j 



= {iy-l)j- 



Using these identities we find 



D 



u,k = e s 



^ (i- sY ^ fr'^Kj)^ ii-sY+^ ' ii-sY+'' 



(D.5) 



(D.6) 



(D.7) 



The first and the last term can be reabsorbed in the sum, which yields the identity (D.l). 
E Recurrence relation for jTl^^fc 

In this appendix we show that the integral J^u,kj defined on the LHS of eq. (6.5), satisfies 
the recurrence relation 

k 

Ju,k = -6"-^e-"-^(i/)fc + ^(j + l)k-j [abX_2,j + l)J.-i,j] ■ (E.l) 

j=0 

To prove this relation we substitute eq. (6.9) for ^/j^ ^ and rewrite the integral as 



Ju,k = -b' 



v-l 



ds [i 



-bs 



I _^ fk 



e \ e s 



E 



k-j ^v-k+j-l^-^ _ 



v)ja s 
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3=0 

X [a(l -s) + {u-k + j- l)s(l -s)-{k- 3)s^] , (E.2) 

where in the second step we have integrated by parts (noting that only the j = k term 
contributes to the surface term) and in the last step we have performed the derivatives. We 
now manipulate the sum over j in eq. (E.2), which we denote by 5. We split the coefficient 
{ly — k + j — 1) of the second term in — 1 and —{k — j) and merge this second part with 
the third term to obtain 

S = Y,('') {j^)ja''-^s''-''+^-^{l - s)^-^'-i [a(l -s) + {iy- l)s(l - s) - {k - j)s] . (E.3) 

Since for the third term the j = k contribution to the sum vanishes we shift (for the third 
term only) the upper limit of the sum to k — 1 and then j — )• j — 1, resulting in 

k 



S = Y, (^)a'=-^s^-'=+^~3(i - s)fc-i[(j.).a + (z. - l),+is] 
j=i / 

E 0) - s)'-^ [{u - l),a l),+^s] , (E.4) 



j=0 

where in the last step we used the identity 



The Pochhammer symbols satisfy 

{u - 1), = {u- 2) /"^ + ^ ={u- 2), + j{u - l),_i , (E.6) 

and repeating this argument j times we find 

j 

{u - 1), = J2U - i + - 2),-. • (E.7) 

Similarly, 

{V - l),+i = {V- l)j{u -l+j) = [y-l){y- 1)^. 1)^. , (E.8) 



-27- 



which after j repetitions yields 

j 

{u - = {u-l) ^{j - i + lUu - . (E.9) 

i=0 

Using the identity 

repeatedly, we find from eq. (E.7) 

n - 1).- = E i^) + - 2).- = Y.{k-z + l\ 2),-. (E.ll) 

and from eq. (E.9) 

Q) - = - 1) E Q) - ^ + - 

= {v-l) j^ik + ') {v - 1),-.. (E.12) 

i=o ^-^ 

We now substitute the identities (E.ll) and (E.12) in eq. (E.4) and find 

k 



j=0 j=0 ^-^ 

= ^^(A; - ^ + l)i E Z a^-'s^-^^^-^l - s)^-^[{v - 2),_,a + [u - l){u - l),_,s] 

i—n o— n \ J / 



i=0 j=0 



(E.13) 



where in the second step we have interchanged the sums over i and j (see the shaded 
triangle in figure 3) and in the last step we have shifted the index j by i. Looking back at 
the eq. (6.9) we see that J'^^k of eq. (E.2) can be written as 



k 







_^.-ig-a-6(^)^ + Y,ik + [abj,^2,k-i + (z^ - l)Ju-l,k^i] , (E.14) 



i=0 

which, after introducing j = k — i, proves the relation (E.l). 







j 



k 



Figure 3. The range of the indices in the double sums of eq. (E.13) is given by the shaded area. 

F Proof of identity (6.17) 

With the definitions 



Tu,n = abSu^2,n + (z^ - l)S,y-i^r. 
k 

n=0 



(F.l) 
(F.2) 



we need to show that 



^,,k + b''-'{u)k + a''~'kl = S,^k- 



(F.3) 



We substitute the definition (6.3) of S^^k in (F.l) and perform some manipulations that 
are explained at the end of the string of equations, 



I/— 3 v—'i—i 

i=0 j=0 

u-2 u-2-i 



{u -3-i)\{u -3 + n- j)! 
(i/ - 3 - 



+ ('^-l)(^-2-^)!(^-2 + n-j)! ^,^^ 



{ly - 2 - i - j)W.j\ 
(^_2-i)!(z.-3 + n-i)! , 



i=0 j=0 
u-2 v-2-i 

U {u - 2 - ^ - jy.{i - ly.jl 

^ g "^2^ ('^-l)(^-2-0!(^-2 + n-j)! ^,^, 

u-2 



{v - 2 - i - 



u-2 



(^_l)!(^_2 + n-j)! ^ ^ 

^ {v-2-m ^ 



u-2-i 

+ E 

j=0 



'v-2-i 



(i/-2-i)!(z.-3 + n-i)! ^,^^.^i 



i=o 



{u - 2 - i - jy.{i - ly.j] 



(^_l)(^_2-i)!(i.-2 + n-j)!^,^^. 



{i^ - 2 - i - jy.iljl 
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V (^-l)!(^-2 + n-j)! ^ 



j=0 

u-2-i 

+ E 

j=0 
u~2 



1=1 



[y - \)(v - 2 - i)\(y - 2 + n - j)\ 



{y-2-i-j)\i\3\ 



h'o? 



(^_l)!(^_2 + n-i)! ^ sr- 



j=0 



i=l 



n-l + i)!^,^,_i_, , (i.-l)(z.-2 + n)!^ 



(.-1)! 



+ 



+ 2^ + (i^- 1 - ^ - j)) ^ ^ : rwT^. 6 a-' 



(i/ - 1 - i - j)!i!j! 



I/-2 



i=l j=0 



(z^ - 1 - i - j)!z!j! 



(F.4) 



In the second step, we have shifted z — )• z — 1 m the first sum. In the third step, we have 
extracted the i = term of the second sum and merged the remaining sums. In the fourth 
step, we have shifted j — )• j — 1 in the first sum in square brackets. In the fifth step, we 
have extracted the j = v — 1 — i term of the first sum and the j = term of the second 
sum in square brackets and merged the remaining sums. In the last step, we have written 
~ j ~ 2)! as (v — 1 — — 1 — jy. in the first sum. Also, we have observed that 

U + ~ ~ 1 ~ ^ ~ j) = (z^ — 1 — i)(z^ — 1 — j) and merged the three terms in square 
brackets into a single sum. 

The result of eq. (F.4) suggests to merge the two sums since the first sum is just the 
i = term in the second sum, up to a missing i = 0, j = v — l term. However, such a merge 
needs to be done with care as the extra term is formally undefined for n = 0. We therefore 
first split u — 1— j as (i/ — 1 + n — j) — n in both sums and substitute T^^n in eq. (F.2) . We 
switch the order of the summations and first perform the sum over n. Considering only 
the terms that depend on n, the sum to be evaluated is 



^(n + l)fc_„[(i/ - 1 + 



n-j'j 



n{v — 2 + n — j) 



n=0 



J](n + l)fc_„(i/ - 1 + 



n-j) 



'^{n)k-n+i{i^ - 2 + n - j)! 



n=0 

■ k k-r 

E-E 

.n=0 n=0. 



n=l 



(n + l)k-ni'^ -l + n-j)\ = iu-l + k-j)\, 



(F.5) 



where in the second step we have used n(n + l)k-n = (f^)fc-n+i and observed that the 
n = term does not contribute to the second sum, and in the third step we have shifted 
n — >• n + 1 in the second sum so that only the n = k term remains. Using this result we 
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obtain 

j=o ^ ' •' i=l 3=0 ^ J) J 

1=0 j=0 ^ J) J 

h h ^-^-^-m3\ (--1)! ' ^ 

where in the second step we have merged the sums and corrected for the i = 0, j = L' — 1 
term, and in the last step we have extended the sum over i up to — 1 and corrected for 
the z = 1/ — 1, J = term. Looking back at eq. (6.3) and using (i^ — 1 + A;)!/(z^ — 1)! = (i^)^ 
this yields 

^u,k = S,,k - kla""-' - {u)kV'-^ , (F.7) 

which completes the proof. 



G Error estimation for the random matrix simulations 

When using reweighting the ensemble average (7.4) can be written as z = x/y, and the 
error on z is given by the usual error propagation formula 

a, = \-z\M + ^-2p,,y^^, (G.l) 
y X y 

where ax = a^/ViV and ay = ay/^/N are the standard errors of x and y, a^ and ay are 
the square roots of the sample variances of x and y, pxy is the correlation coefficient of x 
and y, and is the sample size. 

In the quenched simulation all configurations are independent and the standard errors 
are computed using the total sample size A^. For the phase-quenched and sign-quenched 
simulations the ensembles are generated using a Metropolis algorithm, and the autocor- 
relations in the Markov chain have to be taken into account by modifying the standard 
errors to 

where Tint,x and T\nt,y are the integrated autocorrelation times for x and y [24]. In eq. (G.l) 
the correlation coefficient pxy is computed over the complete sample, and it is assumed 
that Pxy is not affected by different autocorrelation times for x and y. 
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